#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "galois.h"
#include "galois_add.h"
#include <sys/time.h>

#define K 2
#define L 32
#define LK 64
#define Q 4294967296

int fxn_gcd(uint32_t *two, uint32_t s, uint32_t j)
{
	uint32_t *one, *tmp, shift_count, inv;
	int i, high_bit_1, high_bit_2, array_size;

	if (j == 1) array_size = Q+1;
	else if (j == 2) array_size = Q*Q + 1;

	one = (uint32_t *) malloc(sizeof(uint32_t) * array_size);
	tmp = (uint32_t *) malloc(sizeof(uint32_t) * array_size);

	for (i = 0; i < array_size; i++) {
		one[i] = two[i];
		two[i] = 0;
	}
	two[0] = 1;
	two[1] = s;
	two[2] = 1;

	for (high_bit_1 = array_size-1; high_bit_1 >= 0 && one[high_bit_1] == 0; high_bit_1--);
	for (high_bit_2 = array_size-1; high_bit_2 >= 0 && two[high_bit_2] == 0; high_bit_2--);

	while(1) {
		for (i = 0; i < array_size && two[i] == 0; i++);
		if (i == array_size) break;

		while(high_bit_1 >= high_bit_2) {
			for (i = 0; i < array_size; i++) tmp[i] = two[i];

			shift_count = high_bit_1 - high_bit_2;

			inv = galois_composite_divide_k2s(one[high_bit_1], two[high_bit_2], LK, s);

			for (i = 0; i < array_size; i++) tmp[i] = galois_composite_multiply_k2s(inv, tmp[i], LK, s);
			for (i = 0; i <= high_bit_2; i++) one[i+shift_count] ^= tmp[i];
			for (high_bit_1 = array_size-1; high_bit_1 >= 0 && one[high_bit_1] == 0; high_bit_1--);
		}

		for (i = 0; i < array_size; i++) tmp[i] = one[i];
		for (i = 0; i < array_size; i++) one[i] = two[i];
		for (i = 0; i < array_size; i++) two[i] = tmp[i];

		for (high_bit_1 = array_size-1; high_bit_1 >= 0 && one[high_bit_1] == 0; high_bit_1--);
		for (high_bit_2 = array_size-1; high_bit_2 >= 0 && two[high_bit_2] == 0; high_bit_2--);
	}

	for (high_bit_1 = array_size-1; high_bit_1 >= 0 && one[high_bit_1] == 0; high_bit_1--);

	if (high_bit_1 == 0) return 1;
	else return 0;
}

main(int argc, char **argv)
{
	uint64_t s, one;
	int i, j, reducible, array_size;
	uint32_t k, z;
	uint32_t *two;

	one = 1;

	for (s = 2; s < (one << L); s++) { // f = x^2 + sx + 1
		reducible = 1;
		for (i = 1; i <= K/2; i++) { //only one value for k==2, 2 for k == 4
			/* get (x ^ q ^ i - x) % f */
			if (i == 1) array_size = Q + 1;
			else if (i == 2) array_size = Q*Q + 1;
			two = (uint32_t *) malloc(sizeof(uint32_t) * (array_size));

			for (j = 0; j < array_size; j ++) {
				two[j] = 0;
			}

			two[1] = 1;
			two[array_size-1] = 1;

			k = fxn_gcd(two, s, i);
			free (two);

			if (k == 0) { reducible = 0; break; }
		}

		if (reducible) {
			printf("%lu\n", s);

			for (j = 0; j < 100; j++) {
				k = lrand48() % (one << LK);
				z = galois_composite_inverse_k2s(k, LK, s);
				z = galois_composite_inverse_k2s(z, LK, s);

				if (z != k) printf("err with %lu\n", s);
			}
		}
	}
}
